
/************************************************************************/
/*  									*/
/*  autometh.c								*/
/*									*/
/*  Tom Lookabaugh							*/
/*  August 13, 1987							*/
/*  Last Revision: August 13, 1987					*/
/*									*/
/*  DESCRIPTION								*/
/*	LPC analysis by the autocorrelation method.  Uses a set of	*/
/*	autocorrelation coefficients and computes the autocorrelation	*/
/* 	of the prediction filter coefficients and the prediction 	*/
/*	error.  Note that m is the order of the prediction; hence	*/
/*	r has dimension m + 1.  This routine is based on the		*/
/*	Levinson's algorithm implementation of Figure 3.1, p. 51 of	*/
/* 	J. D. Markel and A. H. Gray, Jr., "Linear Prediction of 	*/
/*	Speech," Berlin: Springer-Verlag, 1976.				*/
/*  RETURNS								*/
/*	No significant return values.  Modifies storage pointed to 	*/
/*	by ra and alpha.						*/
/*  CALLS								*/
/*	error(), warning()						*/
/*									*/
/************************************************************************/

#include "lib.h"


void autometh (r, m, ra, alpha) 

    double *r;		/* autocorrelation coefficients	 (size m + 1)	*/
    int m;		/* filter order 				*/
    double *ra;		/* autocorrelation of best filter (size m + 1)	*/
    double *alpha;	/* prediction error				*/
    
{


/* See Markel and Gray for a description of the internal variables.  Note
   that their code has been modified by (i) assuming the autocorrelation
   coefficients were calculated externally (eliminates code upto their 
   label 100), (ii) allocating internal memory for rc, b, and a dynamically,
   and (iii), decrementing pointers to allow for indexing differences 
   between C and FORTRAN.						*/

    double  *b, s, *a, *rc;
    int mf, minc, j, jb, ip, mip, k;



    if (((rc = (double *)calloc(m, sizeof(double))) == NULL) ||
	((b = (double *)calloc(m, sizeof(double))) == NULL) ||
	((a = (double *)calloc(m + 1, sizeof(double))) == NULL))
      error("not enough memory (autometh 1)");
    r--;
    a--;
    ra--;
    rc--;
    b--;



    *alpha = r[1];
    rc[1] = -r[2]/r[1];
    a[1] = 1.0;
    a[2] = rc[1];
    *alpha = *alpha - rc[1]*rc[1]*(*alpha);
    mf = m;
    for (minc = 2; minc <= mf; minc++) {
        m = minc - 1;
        for (j = 1; j <= minc; j++) {
            jb = minc - j + 1;
            b[j] = a[jb];
        }
        m = m + 1;
        s = 0.0;
        for (ip = 1; ip <= m; ip++) {
            mip = m - ip + 2;
            s = s + r[mip]*a[ip];
	}
        rc[m] = -s/(*alpha);
        for (ip = 2; ip <= m; ip++) a[ip] = a[ip] + rc[m]*b[ip - 1];
        a[m + 1] = rc[m];
        *alpha = *alpha - rc[m]*rc[m]*(*alpha);
        if (*alpha <= 0) warning("insufficient accuracy (autometh 2)");
    }

    for (j = 1; j <= m + 1; j++) {
        ra[j] = 0.0;
        for (k = 1; k <= m + 2 - j; k++) ra[j] += a[k]*a[k + j - 1];
    }

    b++; rc++; a++;
    free(b);
    free(a);
    free(rc);
}

